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Abstract 

Classical solvable stochastic volatility models (SVM) use a CEV process for 
instantaneous variance where the CEV parameter 7 takes just few values: 
- the Ornstein-Uhlenbeck process, 1/2 - the Heston (or square root) process, 1- 



GARCH, and 3/2 - the 3/2 model. Some other models were discovered in Henry 



Labordere ( 2009| ) by making connection between stochastic volatility and solvable 



diffusion processes in quantum mechanics. In particular, he used to build a bridge 
between solvable (super) potentials (the Natanzon (super) potentials, which allow 
reduction of a Schrodinger equation to a Gauss confluent hypergeometric equa- 
tion) and existing SVM. In this paper we discuss another approach to extend the 
class of solvable SVM in terms of hypergeometric functions. Thus obtained new 
models could be useful for pricing volatility derivatives (variance and volatility 
swaps, moment swaps). 



*I thank Peter Carr, Mark Craddock, Alexander Antonov, Igor Halperin, Roza Galeeva and Michael 
Spector for useful comments and discussion. I assume full responsibility for any remaining errors. 



1 Introduction 



Analytical tractability of any financial model is an important feature. Existence of 
a closed-form solution definitely helps in pricing financial instruments and calibrating 
the model to market data. It also helps to verify the model assumptions, check its 
asymptotic behavior and explain causality. In mathematical finance many models were 
proposed, first based on their tractability, and only then by making another argument. 

Stochastic volatility models (SVM) were introduced to resolve shortcomings of the 
Black-Scholes model. They are highly used to evaluate derivative securities, such as 
equity and FX options, and variance/volatility products such as variance/volatility 
swaps. Tractability of these models is limited to some partial cases (which for this 
reason became very popular and play the same role for SVM as the Black-Scholes 
model for local volatility models). Classical solvable SVM use a CEV process for the 
instantaneous variance where the CEV parameter 7 takes just few values: - the 
Ornstein-Uhlenbeck process, 1/2 - the Heston (or square root) process, 1- GARCH 
or the geometric model, and 3/2 - the 3/2 model. Closed- form solutions for the op- 
tion prices written on the underlying spot were provided by using Fourier inversion if 
characteristic function E[e mJ5ft ] of the underlying process X t is known in closed form. 
The latter is trivial to find for 7 = 0, and for 7 = 1/2 it was given in Heston (1993). 
Later the cases 7 = 1/2, 1, 3/2 were investigated by Lewis (2000) who for this purpose 
developed a method of the generalized Fourier transform. Using that Lewis derived a 
nice representation for the characteristic functions in the above cases, while at 7 = 1 it 
is quite complicated and expressed via series of the Pochhammer indices. This expres- 
sion could be further simplified if there is no correlation between the spot and variance 
Brownian motions. 

For pricing variance and volatility derivatives such as variance swaps and options 
on them one has to know either a characteristic function of the process, or a Laplace 
transform of the quadratic variation of the process E[ e - A M]. For SVM [X) t = £v t dt, 
where Vt is the instantaneous variance. For instance, under the Heston model a fair 
price of the variance swaps was widely reported in the literature (see, for instance, 



Swishchuk (2004)), while for the 3/2 model it was obtained in Carr and Sun (2007) 



(see also Itkin and Carr (2010)). Carr and Sun also derived a closed form expression 
for the joint Fourier-Laplace transform for the 3/2 model. 

Some other volatility derivatives could be priced if the characteristic function is 
known in closed form, for instance gamma swaps (see Lee (2008)), but not the volatil- 
ity swaps (for those only some approximations are available, see Gatheral (2004), 



Swishchuk (2004) ) and high moment swaps (Schoutens (2005)). Same is true if the 
Laplace transform of the quadratic variation is known (Carr and Sun (2007)). For 
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pricing options on variance it is sufficient to have the characteristic function of the 
underlying process, and then use FFT ( |Lee| ( |2004| >, |Sepp| ( |2008| )). 

The following observation is important for calibrating term structure of, e.g. the 
variance swaps to market data using the above models. Consider the Heston model 
that uses a CIR process for the instantaneous variance which is linear in drift. This 
results in the fair price of the variance swap to be independent of how the volatility of 
variance is specified. In contrast, when the drift of the instantaneous variance process is 
nonlinear (e.g. quadratic as in the 3/2 model) the variance swap price does depend on 
volatility of variance as well as on correlation between the spot and the instantaneous 
variance. This provides better flexibility of the model at calibration, and makes it 
richer in terms of modeling variance swaps term structure. 

In addition to the above mentioned solvable SVM some other models were discov- 
ered later in Henry-Labordere (2009) by making connection between stochastic volatil- 
ity and solvable diffusion processes in quantum mechanics. In particular, he used 
to build a bridge between solvable (super)potentials (the Natanzon (super)potentials, 
which allow reduction of the Schrodinger equation to the Gauss confluent hypergeo- 
metric equation) and the existing SVM. Two-dimensional Kolmogorov equation has 
been converted into the Schrodinger equation with a scalar potential, and it was shown 
that the Heston model, the geometric Brownian model (GARCH) and the 3/2 model 
belong to the Natanzon class. Henry-Labordere (2009) also claims a new result that 
the case 7 = 2 (more rigorously volatility of volatility a (a) = a + 77a 2 , rj = const) can 
be solved in closed form in terms of hypergeometric functions, but to the best of our 
knowledge the result has not been presented. 

Despite a high popularity of these models, in the first place because of their 



tractability, in general they are not able to explain many empirical observations. Carr 



and Sun (2007) give a short survey of the existing literature which provides an em- 
pirical support, that being extracted from the market data, 7 > 1. |Carr and Sun 



(2007) consider two groups of papers: one which deals with statistical processes, and 
the other one which considers statistical and risk-neutral processes. In the first group 
using affine drift Ishida and Engle (2002) estimate the 7 to be 1.71 for S&P500 daily 
returns measured over a 30- year period. Javaheri (2004) also estimates this process on 
the time series of S&P500 daily returns, but with the CEV power constrained to either 
be 0.5, 1.0, or 1.5. He concludes that a power of 1.5 outperforms the other two possible 
choices. Chacko and Viceira ( 1999 ) use spectral GMM estimation on the same process 



as in Ishida and Engle. Using the CRSP value weighted portfolio, they estimate the 
CEV power at 1.10 using weekly data over a 35-year period, and at 1.65 using monthly 
data over a 71-year period. In the second group Jones (2003) examines daily S&P100 
returns and implied volatilities over a 14-year period. Using this data with the statis- 
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tical version of the affine drift CEV process, and for S&P100 daily returns, he finds 
the CEV power to be 1.33. For shorter maturity options he concludes that jumps may 



be needed. Bakshi et al. (2004) look at time series of S&P100 implied volatilities as 



captured by VIX. They find that a linear drift model is rejected in favor of a nonlinear 
drift model, and 7 = 1.27. They conclude that 7 > 1 is needed to match the time 
series properties of the VIX index with the CEV models. 

The above results tell us that 7 is not a constant and varies depending on the 
market data, time to expiration etc. So it would be desirable to have 7 to be a 
parameter of the model. Also a nonlinear drift often better fits the market data than 
its liner counterpart. Finally, it would be feasible to have the drift being a function of 
time to calibrate the model to the existing term-structure of the market data (like in 
short-term interest rates models). 

Unfortunately, we do not believe one can obey all these requirements and have a 
tractable model at the same time. Therefore, we need a compromise. In this paper 
we partly sacrifice by tractability in a sense that the characteristic function or the 
Laplace transform of the quadratic variation is not required to be known in closed form. 
However, marginal density of the instantaneous variance process should be available 
in closed form. If this is the case, after the model is calibrated, prices of all volatility 
derivatives (variance swaps, volatility swaps, higher moment swaps, VIX futures, etc) 
can be computed analytically (or semi-analytically, i.e. in quadratures). 

To calibrate such a model, e.g. to vanilla option prices, two approaches seem to 
be reasonable. The first one deals with the numerical pricing. To solve a correspond- 
ing pricing PDE finite differences together with a very efficient splitting method (e.g. 



In't Hout and Welfert (2007)) could be used. The total complexity of the method is 



O(NxM) N - the number of mesh points in the spot space, M - the number of mesh 
points in the variance space. 

The second approach is almost analytical and makes use of the mixing theorem. 
Conditional on the whole path of the instantaneous variance from vq to Vt the option 
price Ci then is given by the Black-Scholes formula with an efficient spot price S e and 
an efficient volatility a e (see Romano and Touzi ( 1997[ )). This method is very fast 



due its almost analytical nature, and we recommend it for calibration of the proposed 
models. 

The rest of the paper is organized as follows. In the next section we specify a 
general SVM which fits the above requirements and derive a backward Kolmogorov 
PDE for the density of the instantaneous variance. To solve it in Section 3 we use a 
Lie symmetries method to find all Lie symmetries and groups of transformations of the 
Kolmogorov PDE. Using the ideas of Craddock (2009) we show that the density which 
solves the Kolmogorov equation could be represented as an inverse Laplace transform 
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of a known function. In Section 4 we show under which assumptions on time- dependent 
coefficients of the model this inverse Laplace transform could be calculated analytically. 
We also derive closed form solutions for the volatility derivatives by computing the 
corresponding expectations. Next Section describes various approaches to the model 
calibration and presents results of some numerical tests. The last Section concludes. 

The proposed model and closed form solutions for the volatility derivatives are new 
and constitute the main contribution of this paper. 

It is important to note^that the model derived in this paper to describe dynamics of 



the instantaneous variance looks similar to the model which was proposed in Carr and 



Linetsky (2006) to describe dynamics of the underlying spot process for the valuation 



of corporate liabilities, credit derivatives, and equity derivatives. The differences are 
as follows: Carr and Linetsky (2006) use a different method (theory of scaled Bessel 
processes) in order to compute some moments (via expectations) of the underlying 
process. They did not provide a density function in the closed form. The main area of 
applications of their results are credit and equity derivatives. In our paper the main 
area are volatility derivatives and options on them. We use a different mathematical 
approach to find the analytical representation for the density of the instantaneous 
variance process and all moments, including higher moments. As a consequence, we 
are able to price in the closed form all the above mentioned volatility products. 



2 Specification of the SVM 

We want to analyze a rather general form of the SVM which allows pricing of path- 
independent contingent claims such as volatility and variance swaps and options on 
them. Consider a frictionless market for the underlying spot price S t ,t G [0, T] for 
maturity T. Under no arbitrage there exists a risk-neutral measure Q such that the 
prices of all non-dividend paying assets are martingales under this measure. The risk- 
neutral process for the underlying price is: 

dS t = (r- q)S t dt + S t ^T t dW } t e [0, T), 

where dW is a Q standard Brownian motion, vt is the instantaneous variance, r, q are 
the interest and continuous dividend rates. Assume that the risk-neutral process for vt 
is given by 

dv t = [q{t)Q{v t ) + s{t)S{v t )} dt + G{t,v t )dZ t , t e [0,T], 

J This was kindly noticed by P. Carr 
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where Z t is a Q standard Brownian motion, whose increments have known constant 
correlation p G [—1, 1] with increments in the Q standard Brownian motion Wt, i.e.: 



E Q [dZ t ,dW t ]= pdt, te[0,T]. 

Here q(t),s(t) are some deterministic functions of time, which assumed to be a 
model input, Q(v t ), S(v t ) are some yet unknown functions of v t . As it was discussed 
earlier, using q(t), s(t) to be a function of time, versus to be a constant as in a standard 
specification of the SVM, allows a better calibration of the model. This is especially 
important if the model is used for calibration of the term structure of the variance 
swaps. 

Note, that this is rather general representation of the SVM. Let us outline which 
particular functional form of q(t),Q(v t ), s(t), S(v t ),G(t,v t ) admits a closed-form solu- 
tion that was discussed in Introduction. To the best of our knowledge only a CEV 
process has been used in the literature for G(t,v t ), i.e. G(t,v t ) = t> 7 , 7 > 0. Under 
this model the following solutions are known. 



1. 7 = 1/2, q(t) = K6(t),Q(v t ) = l,s(t) = K,S(v t ) = v t - the Heston model (Heston 



(1993)). 



2. 7 = l,q(t) = K9(t),Q(v t ) = l,s(t) = K,S(v t ) = v t - GARCH-type model. 



Lewis (2000) derived the characteristic function for this model which is rather 



complicated. It can be simplified, if p = 0. 



7 = 3/2, q(t) = K6(t),Q(v t ) = l,s(t) = K,S(v t ) = v t - 3/2 model. Lew5|(2000) 



derived the characteristic function for this model. Later Carr and Sun (2007) 
extended this case and obtained the closed form solution for the characteristic 
function and Laplace transform of the quadratic variation when Q(vt) = v t , s(t) = 

8,S(v t ) = vl 



4. 7 = 2,q(t) = K6(t),Q(v t ) = 1, s(t) = n,S(v t ) = v t - |Henry-Labordere| fl2009| ) 
claims that this model can be solved in closed form in terms of hypergeometric 
functions, but the result has not been presented. 

Let us also remind those features of the model that are desirable in our setup. 

• We still want to use a CEV model of volatility of volatility, e.g. G(t, v t ) oc v 7 , 7 > 
0, but with 7 being a calibrated model parameter. 

• Specification of q(t), s(t), Q(v t ), S(v t ) should preserve mean-reversion. 
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• To have the model being suitable for calibration of the term structure we want 
to preserve time dependence of, at least, some coefficients, e.g. q(t) or s(t), or 
better both. 

Based on these requirements to the proposed model consider the following SDE for 
the instantaneous volatility Xt = ^fvt 



(1) 



dx t = [q{t)x a t - s(t)xf\ dt + l(t)x] +1 dZ t 
Application of the Ito's lemma gives the corresponding SDE for v t 



dv t 



2q{t)v t 2 + l 2 {t)v2 +l - 2s(t)v t 2 dt + 2l(t)v t 2 dZ, 



6+1 



7+2 



(2) 



a ^ b, are some values to be determined, l(t) is some deterministic 

7+2 



Here a, b G 

function of t. When deriving the Eq. (2) we assumed that G(t,v t ) = 2l{t)v t ' , 7 > 
0, l(t) > 0. Based on the sign of s(t) and q(t) it could be that a > b or b > a to 
preserve mean reversion. 

Our further goal is to find under which values of a and b, and particular functional 
form of q(t),s(t) and l(t) the risk-neutral density of the v t could be found in closed 
form. This is not a trivial problem as it could seem at the first glance. Indeed, for a 
standard CEV process 

dSt = fxStdt + 5Sj +1 dZ t 

with the drift fi = const and parameters 5, 7 = const it is known that change of variable 
1/(<5|7|)S7 7 reduces the CEV process without drift (// = 0) to the standard Bessel 



z t 



process of order I/27 (see Davydov and Linetsky (2001), Revuz and Yor ( 1999[ )). Then 
the continuous part of the risk-neutral density of St, conditional on So = S, is obtained 
from the well known expression for transition density of the Bessel process. If \x 7^ 0, 



using the result of Goldenberg (1991) this CEV process could be obtained from the 



process without drift via a scale and time change 



S$ = c^S. 



T(*)» 



T{t) 



1 



2/i 7 



3 2jt7t 



1). 



Carr and Linetsky (2006) further extended this approach. They model the price of 



the defaultable stock under an equivalent martingale measure as a time-inhomogeneous 
diffusion process S^ , t > with state space E A = (0, oo)Z7A. Here if the process S t hits 
zero it is sent to the cemetery state A at the first hitting time of zero, T . Specification 
of the model looks similar to the Eq. ([T]), namely: if x is the initial value Sq = x > 0, 
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then the diffusion coefficient is a(t)x 7+1 , the drift is [r(t) — q(t) + b(t) + ca 2 {t)x 2 ' 1 ]x, 
where r(t),q(t),b(t),a(t) are some functions of t. Carr and Linetsky (2006) call this 



stock price process as the jump to default extended CEV process, or JDCEV. They 
also give a survey of predecessor's papers that considered a similar model with constant 
coefficients a, b, c. 

Using a theory of scaled Bessel processes they further managed to derive closed 
form valuation formulas for corporate bonds, credit default swaps, stock options, and 
other credit and equity derivatives. 

In this paper we, however, use a different approach. Main reason is that to fulfil our 
program as it was stated above, we need to find the density of the underlying process 



(instantaneous volatility x t ) in closed form, that is not provided in Carr and Linetsky 
( |2006| ). 

To finalize description of the model it is important to note that for 7 > according 
to Feller's classification the origin vt = is a natural boundary, and infinity is an 
entrance boundary 



sec 



Davydov and Linetsky (2001)). 



3 Backward equation for the density and Lie sym- 
metries 

Under the SVM with no jumps various volatility derivatives can be represented as 
expectations under the risk neutral measure Q. Then, for instance, fair strike of the 
variance swap (or annualized total expected realized variance) conditional on the initial 
level v = [v t \t = 0] reads 



V 2 = E< 



T 



v s ds \vq = v 



1 

T 



EqK \vq = v]ds. 



Similar expression gives a fair volatility swap value 



1 



T 



T 



v s ds \vq = v 



(3) 



European call option price written on the underlying realized variance (e.g. VIX 
options) is by definition 



Co 



-rT 



T 



x s ds 



K 



X = X 



S 



with the standard notation (y — K) + = max(y — K, 0). 
For the put option a similar expression reads 



-rT 



1 

T 



x Q ds 



X 



Let p(x, t : i', y) be a backward transition density from the state y at time t' to the 
state x at time t. Let v = t' — t be a backward time. A standard argument tells us 
that the expectation u(x, v) = E[$(x) \xf — x] solves the following Cauchy problem 

dU ^ U) = [*(„)*■ - s(u)x b ] 9 -^l + \vfx^ d2 <Y\ (4) 
av L J ox 2 ax 2 

0) = $(x t >) 

This PDE is a backward Kolmogorov equation. The transition density p(x,t : t',y) 
is also a fundamental solution of the Eq. Q. Thus, if we know the density, the 
expectation could be computed according to u(x,t') = J °° §(y)p(x, : t',y)dy where 
it was assumed that we price all products at t — 0. 

However, the fundamental solutions are not unique. In contrast, the probability 
transition density p(x,t : t',y) for the process x t is unique and obeys the additional 
condition J °°p(x, t : t', y)dy = 1. Therefore, it is not efficient to use standard methods 
to solve the Eq. Q if one needs to determine the density. Indeed, suppose we use a 
change of variables method that reduces the Eq. ([5]) to, say, the heat equation. The 
density for the heat equation is known. However, under backward transformation to 
the original variables, this function will not integrate to 1, despite it still remains to 
be the fundamental solution of the Eq. Q. Therefore, it is difficult to distinguish the 
density from the other solutions. Moreover, as shown in Craddock and Lennox (2007) 



this requires a theory of generalized functions and distributions. 



Recently Craddock (2009) proposed a new method to find the density, which utilizes 



Lie symmetry analysis of the parabolic PDEs. What actually Craddock showed in his 
paper is that for the PDEs with nontrivial Lie symmetry algebras, the Lie symmetries 
naturally yield Fourier and Laplace transforms of fundamental solutions. Therefore, 
our further goal is to derive an explicit representation of such transforms in terms of 
the coefficients of the PDE, and find the solutions that obey the desirable properties 
of our model. 



3.1 Lie symmetries 

A symmetry of a differential equation is a transformation which maps solutions to 
solutions. In the 1880s Lie developed a technique for systematically determining all 
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groups of point symmetries for systems of differential equations. By Lies method, we 
look for infinitesimal symmetries of the form (see 



v = f (x, t, u)<9 x + r(x, t, u)<9 t + 0(x, t, u)d u . 



(5) 



We need to find conditions on £, r, which guarantee that v generates a symmetry of 
the Eq. Q. Standard arguments show that r can only depend on t, and £ can only 
depend on x and £. Further, <p must be linear in it, e.g. <f)(x,t,u) = a(x,t)u. Lie's 
Theorem says that v generates a local group of symmetries if and only if 



pr 2 v 



du 
dt 



[q(t)x a - s(t)x b ] 



du 

<9x 



1 l2 m 2( 7 +l) ^ U 

2 1 ltjX «9x 2 



0. 



(6) 



where pr 2 w is the second prolongation of v, and u is the solution of the Eq. Q. The 
explicit prolongation formula for a vector field is given in Olver (1993). 
If 7 7^ 1 this leads to 



4>(x, t) = a(x, t)u + (3(x, t), 



(7) 



i(x,t)=x^H(t) \a(t) + 



x-i [2r{t)l'{t) + l{t)T'{t)] 
27/(t) 2 



where a(x, t), /3(x, t) are some solutions of the our backward Kolmogorov equation. 
Choose (3(x,t) = 0. Then the only nontrivial solution for a(x,t) could be obtained if 
a = 1, b = 27 + 1, a(x,t) is given by the following formula 

a (x, t) = a + v(t> m , v(t) = c 2 e / ° m?(y) ^, 

where m is a parameter to be chosen, and Ci, c 2 are the integration constants, and s(t) 
reads 

s{t)= l -{m-l)l{tf. 
Substitution of these values into the Eq. ([!]) transforms it to 

1 



dxt 



q(t)x t — -(m — l)l(t) x t 



2^,27+1 



dt + l(t)xl +1 dZ-, 



Now a tedious algebra (that we omit here) shows that the Eq. ^ could be solved 
in 3 cases: m = —27, —7,7. Accordingly, all basis triplets t, u), r(x,t,u), a(x,t,u) 
that produce a vector field in the Eq. ^ and solve the Eq. ^ are given below 

2 In this section for a better readability we revert our notation for the backward time v back to t 
since this should not create any problem. 
3 Remember that <p(x,t,u) = a(x,t)u 
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3.1.1 m = —7 

In this case the explicit form of our model follows from the Eq. (|8 



dx t 



l -(l + l)l{tfx 2 r l + q{t) 



x t 



dt + l(t)xJ +1 dZ t . 



(9) 



Lie algebra of infinitesimal symmetries of the Eq. Q is now spanned by the vector 
fields 



vi = ud u 

v 2 = <yZ(t) n 

v 3 = Z\t) 



i(yf 



dy] x 1+7 cL + =—x x 7 w<9 u 



(10) 



i( t ) 



[d t - q{t)xd x 



v 4 = —x 



27 



q{t)K{t) 



Z(t) 



l(t) 



i(y) 



Z(y) 



dy 



d x + K(t)d t , K(t) = 
v 5 = Z(t)l{0)x 1+ ^d x 

Here Z(t) = e^o ~fq(y)dy_ Among these solutions we are interested just in non-trivial ones 
with a(x,t) 7^ and £(x,t) 7^ 0, i.e. those in line 2 of the Eq. (10). 

According to the Lie's method (see Chapter 2 of Olver (1993)) given the symmetry 
Vi in the form of the Eq. ^ we need to exponentiate it in order to find the one- 
parameter group Gi generated by the Vi. 

The first symmetry in the Eq. (10) is trivial and says that the new solution of the 
Eq. Q can be obtained from another solution by multiplying it by a constant. The 
symmetries in lines 3,4 are hard to reproduce in closed form, since exponentiation can 
not be done for the arbitrary time-dependent functions q(t) and l(t]n The symmetry 
in line 5 translates just x, so 



: {[x^ + lf iZ(t)l(0)\ 



-1/7 



,t,u 



Exponentiating line 2 of the Eq. (10) we obtain 



G 2 : { [x^ + 1 fiUt)]' 1/ \t J exp 



Z(t) 



u 



'At the best the result can be expressed via inverse functions 
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where £ c (t) = ^Z(t) y J * Z(y) 2 l(y) 2 dyj , u is some solution of the Eq. 
group parameter. 

Existence of this symmetry group implies that if u = G(x, t) is a solution of the 
Eq. Q so are the functions 

u»{x,t) = exp j-x-^Z(i)- 1 + ]pZ{t)-%{t)A G {[x- 1 + 1 ^ c (t)Y 1 '\t) ■ 

Note, that the corresponding SDE for the instantaneous variance can be derived from 
the Eq. g using Ito's lemma 

7+2 

dv t = [(2 + 7 )/(t)\ 1+7 + 2q{t)v t ] dt + 2l(t)v t 2 dZ t . 

The first term in the drift is positive for if 7 > —2, therefore to preserve mean reversion 
we must have q(t) < 0. Also v] +1 grows faster than Vf if 7 > 0, therefore the ratio 
q(t)/l(t) 2 has to be a rapidly increasing function of time to compensate this effect and 
provide a mean reverting behavior of the drift. Finally, based on our earlier discussion 
of the empirical data we expect the CEV exponent to vary from 1 to 2, which means 
< 7 < 2. 



4 ) , and \i is the 



Specific form of q(t) and l{t). Despite we are not able to analytically exponentiate 
all symmetries found in above, we can extend the class of tractable symmetries by 
choosing some particular form of q(t) or l(t). Assume that 

1 ('(«) 



9(t) - 7 m ■ 

Under this assumption our models takes the form 



dx f 



-x 



1+27 



(l+ 7 )/W 2 + 



xl'(t) 



dt + l(t)xl +l dZ t . 



(12) 



(13) 



The first term of the drift is positive. Therefore, to preserve mean reversion we 
need l(t) to obey the condition V(t)/l(t) = —6H(t) where 9 is a calibration constant, 
and H(t) is some increasing function of t. To obey this condition we choose 



l(t) 



ee 



- Jo 6H{y)dy 



(14) 



where e is a constant volatility of volatility. With this expression for l(t) our model for 
Vt now reads 



dv t 



(2 + 7)ZWV +7 " 2 



6H(t) 

7 



-v t 



dt + 2l(t)vt* 2 dZ t 



12 



Again two concurrent effects affect the drift. On the one hand v 1+1 grows faster 
than v if 7 > that does not allow mean-reversion. On the other hand the exponent 
in the first term makes it rapidly decreasing if the function H(t) grows fast with time. 
Thus, one can always properly choose H(t) to guarantee mean reversion. 



With functions q(t),l(t) which obey the Eq. (12), (14) symmetry v 4 in the Eq. (10) 
reads 

26H(t)t - 1 
V4 — - 



2 7 



-xd x + td t 



that now could be exponentiated in closed form. 

Some other choices of the function l(t) are also possible for this purpose, for instance 



that gives v 4 



-x 



l(t) 



h + 



exp 



1 

2t 



t 

- I q(z)dz 
'0 



2l 



d x + t 2 d t , etc. 



3.1.2 m = 7 

In this case the explicit form of our model is 

"1 



dx f 



: (l- 7 )/(t) 2 ^ +1 + g(t)a; t 



and 



dt + l(t)x1 +1 dZ u 



7+2 



dv t = [(2 - -i)l(tfv] +1 + 2q(t)vt] dt + 2l(t)v t 2 dZ, 



(15) 



Therefore, again to preserve mean-reversion either q(t) must be negative and 7 < 2 
whereas l(t)/q(t) must be a rapidly decreasing function of t, or 7 > 2 and q(t) > 0. 

Lie algebra of infinitesimal symmetries of the Eq. ^ in this case is spanned by the 
vector fields 



Vi 

v 2 
v 3 
v 4 



ud u 

m 

7 

—x 



x 1+ W x + Z{t)x^ud u 
1 



2 7 



+ q(t)K(t) 



d x + n(t)d t 



-Z(tf 



l(t) 



q{t)xd x + n(t)d t 
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Exponentiating line 2 we obtain 

G 2 : | [x^ + nZ(t)] ~ lh , t, [1 + fiZ(t)x*] u 



(16) 



3.1.3 m = -27 

In this case the explicit form of our model is 

1 



dxt 



il + 2 1 )l{tfxV +1 + q{t)x t 



dt + l(t)xJ +1 dZ, 



(17) 



and 



7 + 2 



dv t = 2 [(1 + 7 )Z(t) \ 1+7 + q{t)v t ] dt + 2l{t)v t 2 dZ t . (18) 

Lie algebra of infinitesimal symmetries of the Eq. Q here is spanned by the vector 
fields 



Vi = ud u 



v 2 = 27 



Kl(*) 



Kv) 
z(y) 



dy + 2yq(t)K 1 (t) 



v 3 



-x 



l(t) 2 

1 



xd x - 4 7 2 /t 1 (t)<9 t + Z(t)- 2 a;~ 27 u9 u 
l(u) 



Z(u) 



du dz 



27 



+ ?(*)«(*) 



<9 X + «(£)<9 t 



v 4 = Z(t) 5 



Z(0) 

L*(*)J 



-q{t)xd x + <9 t 



Only line 2 has a particular interest for our purposes (see next section) since it 
provides a non-trivial transformation of the solution, not just the coordinates. Unfor- 
tunately, it can not be exponentiated at arbitrary l(t). and q(t). Therefore, we will use 



their specific form proposed in the Eq. (12) and (14). Then the vector field V2 reduces 
to 

v 2 = -2t 7 e 2 [-l + t6H(t)}xd x - 2y 2 e 2 t 2 d t + e 2 ^ m{ - y)dy x~^ud u 

To exponentiate it, choose the simplest model H(t) = 7. Thus, in this model the 
speed of mean reversion is determined by the only parameter 9. Exponentiating, we 
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find some new solutions for our problem 



= exp 



1 + 2t 7 2 /i 



(19) 



• G ( e K-^), (1 + 2 , 7V) v^_l_) 



4 Generalized Laplace transform and solutions for 
the density 



The main idea of |Craddock (2009 ) is that as the density p(x, t : y, t') is also the Green's 
function of the backward Kolmogorov PDE, the solution of this PDE can be represented 
in the form 

poo 

Up{x,u)= Ufj,(y,0)p(x,t : y,t')dy (20) 
Jo 

Now let us use the solutions map of the Kolmogorov PDE found in the previous 
section. Few useful observations could be made immediately. 

1. G(x, t) = 1 is the solution of the Kolmogorov PDE given in the Eq. Q. Therefore, 
we plug in this solution into all the symmetry maps from the previous section. 

2. By definition Z{y) = e^o 7<?(y)<fy_ Switching back to time t' and t translates Z(t, t') 
to 

Z(t,t')=e^'^ {y)dy , t'>t. 

Also by obvious reasons only positive values of the initial level of the instanta- 
neous volatility x and the volatility of volatility speed l(t) are considered here. 
Therefore, rewriting £, c {y) by switching back to time t and t' shows that £ c (£, £') 
is positive: 

Ut, t') = 2 1 Z(t ) t') J Z(y, t')- 2 Ky) 2 dy > 0. 

3. At v = (or t = tf) Z(t',tf) = 1, £c{t?,t') = 0. 



'In this section for the backward time t — t' we again use the notation v 
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Let us start with the solutions map found for m = —7 and given in the Eq. (11), 



i.e. with the model given in the Eq. (J9|. Substituting this map into the Eq. (20) and 
taking into account the above observations we obtain. 



-v-^p(x,t:y,t')dy 



-">^Z- 1 (t,t')+w(t,t')ix 2 



(21) 



where w(t,t>) = \^ Z{ y yH{yfdy) > 0. 

But the lhs of this equation is just a generalized Laplace transform of the density, 



exactly in accordance with the main idea of the method in Craddock (2009). 



Few useful theorems justifying this fact are proven in Craddock (2009), where we 



refer the reader to if she is interesting in more details on this method. Below a short 
summary of the results necessary for our further steps is presented. First, to prove that 



the rhs of the Eq. (21 ) is a generalized Laplace transform, observe that this is correct if 
the rhs is a Laplace transform. This follows from a simple change of variables y~ 7 = z. 



Also U/j,(x, t) = e 



-ynZ(t,t')+w(t,t')n 2 



is analytic in and any function analytic in 1/fi 



is automatically a Laplace transform. Thus, m m (x, t) is a generalized Laplace transform 
of some distribution. 

Therefore, we can find fundamental solutions by inverting this generalized Laplace 
transform. Thus found fundamental solution is the transition density. To prove, let 
/i = in the Eq. (21) which gives p(x,t : y,t')dy = 1. 

Now to invert the transform, make change of variables z = y" 1 and rewrite the 



Eq. (21) as 



1 

7 Jo 



-Zfl 



7 + 1 



z ~i p(x, t : z 1 ^ 7 , t')dz 



-X-1 flZ(t,t')+w(t,t')fM 2 



Therefore, 



p(x,t: y,t' 



-( 7 +l)£ 



-x-i nZ(t,t')+w(t,t')ii 2 



(22) 



The Bromwich integral in the rhs is well defined since the integrand vanishes at 
fj, — > ioo and \i — > —ioo. Unfortunately, to the best of our knowledge the inverse 
Laplace transform in the rhs has no closed form representation. 

We can run same machinery for the symmetry map found at m — —27 and presented 



in the Eq. (19), i.e. for the model given in the Eq. (17). This gives rise to the following 



representation of the transition density 



-27 



»p(x,t : y,t')dy = exp 



e 2(V-i)70 x -27„ 

+ 2(t' - t)~i 2 e 2 n 



(23) 
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Inverting this Laplace transform we obtain 



P(x,t : y,t') 



2 7y -(27+l) exp 
e (t' -t)-y9 x --y y-y 



Xt-t)i0 x -*y _j_ -2 7 



_2(t'-i) 7 2 e 2 



2(f - t) 7 2 e 2 

g(*'-*)7^-7y~7 



(tf - t) 7 2 e 2 



(24) 



where is the modified Bessel function of order d. 

To verify that this is the density let us check that p(x, t : y, £') given by this 
expression integrates to 1. Omitting an intermediate algebra we find 



p(x, t : y, t')dy = 1 — exp 



e 2{t'-t)j8 x -2-y 

2{t> - t) 7 2 e 2 



*(0) 



where ^>(x) is the Heaviside theta function, and it is assumed that ^(O) = since we 
work in a space of left-continuous functions. 

On the other hand the direct substitution of \x = into the Eq. (23) verifies that 



p(x,t : y,t') integrates to 1 unless 7 > and x 7^ 0. Therefore, strictly speaking the 
above solution for the density can not be used for 7 > 0. In the next section we will 
discuss this in more detail. 

For the third model found at m = 7 and given in the Eq. (15) we can use the 



solution map in the Eq. (16). Substitution of this map into the Eq. (20) gives rise to 



the following integral equation for the transition density 



y 1 p{x,t : y } t')dy = Z^x 1 



The lhs of this equation is not the Laplace transform. Therefore, we can not use the 
Craddock method to find thus defined transition density. However, one can treat this 
expression either as the Mellin transform^ or just directly guess that the density which 
solves this equation is 



p(x,t : y,t') 



y _ xe ft 



(25) 



3 This was pointed out by A. Antonov 
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5 Pricing variance and volatility swaps 



5.1 Specific form of q(t) and l(t) given in the model Eq. (17) 

(m = —27) 

Since the density of the volatility process xt is known we can price various volatility 
moments of the order i = 1, 2... by computing the corresponding expectation and then, 
if necessary, integrating the result in time. We can do it explicitly in the special cases 
discussed in Section 15. 1.31 

However, as it was mentioned already there are two problems with this solution for 
the transition density. The first one is that at x — > and fi — > the behavior of the 
product x~ 2j fi is undefined if 7 > 0. To highlight the second problem let us formally 



assume that the first problem is resolved, and the Eq. (24) defines the correct transition 



density at 7 > 0. Then the fair price of the z-th moment reads 



Vi 



1 

T 



Q(0,t',x)dt', 



(26) 



Q(t,t',x) 



y l p(x,t : y,t')dy 



2 2 t exp 



r 



2(t' -t)~te x -2-y 



2(f - t) 7 2 e 2 
1 - i + 47' 



2 7 



M 



+ 2{t' - 1) 7 

1 - i + 47 

27 ' 



Gcj) L - K T(2- k)M(k,2,-G), 



G 



W-Wx-^/t, ,<P = 2(t'-t) 1 2 e 2 



[(? - t) 7 2 e 2 ] ^ 



e 2(t'-t) 7 6» a ,-2 7 

2(t' - t) 7 2 e 2 



x 



-2", 



i - 1 



Here T(x) is the gamma function, M(a,b,x) is the Kummer confluent hypergeometric 
function. 

Therefore, in the proposed SVM the price, e.g. of the variance swap, can be ex- 
pressed via hypergeometric functions, similar to what was discussed in Introduction, 



and in a more general case in Albanese et al. (2001). Also as shown in Itkin and Carr 



(2010) under the 3/2 model the variance swaps price has the closed form representation 
also in terms of the confluent hypergeometric function. 

The internal integral Q(0,t',x) exists under some conditions on i. In particular, 
when 7 > it always exists at i — 1. If i — 2 the existence condition requires 
7 > 1/2. Atz = 3 the condition is 7 > 1 and at i — 4 the condition reads 7 > 3/2. 
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i 


Range 


1 


(3 > 1 


2 


(3 > 1.25 


3 


/3 > 1.5 


4 


/3 > 1.75 



Table 1: Ranges of the CEV exponent f3 = (2 + r y)/2 where the Eq. (26) can be used 
to compute the i-th moment of the instantaneous volatility pdf. 



Therefore, these conditions being translated to the CEV exponent (3 = (2 + j)/2 in 
the Eq. (18) mean that first four moments of the pdf for the instantaneous volatility 
can be computed using the above expression, if (3 varies in the ranges given in Tab [TJ 

Also note, that the expression for Vi in the Eq. (26) is well-behaved at t — > t'. 
Expanding the gamma and Kummer functions into series at t — > t' it is possible to 
show that in this limit Vi — > x l ~ 2l ~ l . 

The problem with this solution is that y{G) = GM(k, 2, —G) is an increasing 
function of G since k < 0. As G oc x~ 21 this means that Vi is the decreasing function 
of x. This fact is in contradiction with a usual observation that the variance and 
volatility swaps price increases with the increase of the initial level x. 

This problem can be resolved, however, if we assume that 7 is negative. Based on 
the Eq. (18) suitable values of 7 lies in the range — 1 < 7 < 0. This has few advantages. 
First we can relax the conditions issued on the explicit form of l(t) because with these 
values of 7 the variance equation Eq. (18) is always mean-reverting as far as q(t) < 0. 
Second, this range of 7 eliminates any restrictions on the existence of the solution in 
the Eq. (26) for the arbitrary i. Third, the price now increases with the increase of the 
initial level x, thus following the observable behavior of these instruments. 

It could seem, however, that a drawback of negative gammas is that under these 
conditions the CEV exponent varies from 1/2 to 1, which is not a favorable region (see 
the discussion in Introduction). On the other hand, we can not completely rely on 
the results in the cited papers for the following reason. In our model the volatility of 
volatility in the Eq. Q depends on both 7 and l(t). This makes some problems when 
trying to find their values by calibration. Indeed, various pairs of 7 and l(t) can produce 
same values of the product l{t)x 1+1 . This makes the calibration to be ambiguous if 
we calibrate the model just to the plain vanilla option. A possible resolution of that 
is to calibrate the model to both plain vanilla options and to variance swaps written 
on the same underlying (e.g. S&P 500 and VIX). This would allow an unambiguous 
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calibration and separation of the effect of l(t) from that of 7 on the price of chosen 



calibration instruments. When computing prices of the variance swaps the Eq. (26) 
can be utilized. The results of such a calibration will be presented elsewhere. 



We can also use the Eq. (17) but applying it not to the instantaneous volatility, 
but to the instantaneous variance. This will allow the CEV exponent to vary from 
to 1, thus extending the previous range. 

To calculate volatility swaps we need to know the expectation in the Eq. ^ which 
can not be computed exactly given just the moments of the instantaneous volatility 



pdf. However, it could be found using approximation from Brockhaus and Long (2000) 



E Q [W] «V^M~ 8] ^ff /2 , V=±j\(0,t',x)dt' 



Since now we want to use the Eq. (17) for the instantaneous variance, this formula 
translates to 

Var[V] 



8Vf /2 



where V l is given in the Eq. Q. As Var[V] = E Q [V 2 } - (E Q [V}) 2 = E Q [V 2 } - V? the 
remaining goal is to find Eq[V^J. This is 



^q[V 2 } = t^J o J Q(0,s,x)Q(0,t,x)dsdt. 

In Fig [Tp] the variance swap and volatility swap prices computed using the above 
approach are given as a function of 7 and x. Here we used a model of l(t) given by the 
Eq. ( [HI ) w ^h H(t) = 7 and 9 = 0.1. Also in these tests we used e = 0.1 and T = 0.5. 
The prices are presented in volatility point divided by 100. 

In Fig [3] same results are given in coordinates (x,T) at 7=-0.6, 9 = 0.3, e = 0.1. 
It is seen that the variance swap price slightly decreases with the increase of T, that 
corresponds to the behavior of many SVM, where the long term run is a deterministic 
function of time, for instance, the Heston model. 



5.2 Solutions for 7 > 0, model in the Eq. (15) (m = 7) 



Since for this model the transition density found in the Eq. (25 ) is just a delta function, 
all volatility moments could be priced in closed form, namely 



Vi 



x- 



T 



.1 

'Jo i(.y) d ydt' 
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Figure 1: Price of the variance swap as the function of 7 and the initial level x at 
T = 0.5 
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Figure 3: Price of the variance swap as the function of T and the initial level x at 7 = 
-0.6. 



Therefore, in this model the price of these derivatives does not depend on vol-of-vol 
l(t) as well as on 7. At given T the price Vj is just a linear function of the i-th power of 
the initial level x. This is partly similar to the SVM with a linear mean-reverting drift 
for the instantaneous variance (like the Heston model) where the VIX price is linear in 
the initial level Vq Lin (2007). 



5.3 General case of arbitrary q(t) and l(t), model in the Eq. 

(m = —7). 

The transition density for the model Eq. ^ was found as an inverse Laplace transform, 
which however does not have a closed form representation. Let us formally substitute 



this density given in the Eq. (22) into the definition of the volatility moments 



V; 



1 

T 



-x-t tJ,Z(0,t')+w(0,t')[i 2 



z->y~ 



dt' 



Since we don't know the inverse Laplace transform under the integral in closed form 
it could seem that this double integral could be computed just numerically. However, 
there is a trick which helps to find a closed form representation of the internal integral. 



The idea of the trick consists in differentiating both parts of the Eq. (21) by \i and 
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then putting ji = 0. If we do this just once the result reads 



y y p(x,t : y,t')dy = x ^Z(t,t f ) 



In other words the lhs of this expression is the expected value Eq[x 7 |xq 



x\ 



Since we need to compute the expectation Eq[x j |x = x], i G Z the above result does 
not seem to be directly translated into the necessary form. Nevertheless, this could be 
done using a theory of fractional derivatives. 



Indeed, following Lavoie et al. (1976) consider Weyl fractional derivative which is 
defined as 

D«_^F(x) = L ; / F(t)(t-x)- a ^dt, (27) 



r( 



-a] 



The beauty of this representation consists in the fact that, e.g. for the exponential 
function differentiation keeps its standard rules regardless whether the derivative order 
is integer or real, i.e. 



D«e b * = b a e bx , Re(b) < 0. 



Going back to the Eq. (21 ) the only term in the lhs which depends on fi is e~ y M . 
Or using the notation of the previous equation, this is e bfl , b = — < 0. Taking the 
fractional derivative of the order a of this expression we expect to obtain cy l e bfl , where 
c is some constant. That gives b a = cy % or a = — 2/7, c = (— 1)~ 1 / 7 . 



Further to compute this fractional derivative of the rhs of the Eq. (21) let us use 
the definition of _ oc in the Eq. (27). Also after taking derivatives let fi = 0. By 



doing so instead of the Eq. (21) we obtain 



1)-^ / y l p(x,t:y,t')dy- ' ' 







r(</7) Jo 



-x-<yZ(t,t')+w(t,t')y 2 yi/j-l 



l dy 



The integral in the rhs part does not exist because w(t,t') > 0. This tells us that 
this particular case should be taken out of consideration since the model is not able to 
produce suitable values of volatility derivatives. 



6 Pricing European options on moment swaps 

European options on various moments of the volatility can be directly calculated if the 
density of the underlying process is known in closed form. Indeed, the price C of, say 
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e~ rT E 



the call option on by definition is 

C n (x,T) 

which can be rewritten as 

e rT C n (x, T) = — 



x n t dt - K 



Xq = x 



K l/r, 



{y 11 - K)p(x,t,y)dy 



dt 



This double integral in general can not be computed in closed form, because of non-zero 
low limit of integration. Therefore, numerical integration is required. However, it can 
be done very efficiently, and also in parallel at any parallel architecture. 

In the particular case of the model Eq. (15) (m = 7), this could be done in closed 
form since the density function is just a delta function as it is presented in the Eq. (25). 
Substituting this into the above formula and integrating we obtain 



e rT C n (x, T) = ax n — K = V n — K 
If, 



(28) 



a 



It is important to underline that European options on volatility swaps can not be 
represented even in integral form under the proposed models and have to be computed 
numerically. 



7 Conclusions 

The main result of this paper is a set of new stochastic volatility models which consider 
the instantaneous variance (volatiltiy) process to be a mean-reverting CEV process. 
However, we do not fix the CEV power 7, but rather treat it as a model parameter to be 
determined by calibration. This approach meets the empirical observations presented 
in the literature and gives more flexibility when calibrating the model parameters to 
the market data on plain vanilla options and volatility derivatives. 

The model preserves mean reversion and time dependence of the mean-reversion 
coefficient which also helps in calibrating term structure of variance swaps and other 
volatility derivatives. It also contains a non-linear drift term. 

What distinguishes our proposed set of models from many possible models of this 
type is that prices of the variance and moment swaps are obtained in closed form 
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(by computation of one integral in time). They also give an approximate solution 
for the volatility swaps. Also options on moments of the realized volatility are priced 
via computation of a double integral. This could be efficiently done using parallel 
calculations. 

Despite our models do not allow a closed form solution for the characteristic or 
moment generations function, still they could be efficiently calibrated, e.g. to vanilla 
European option prices. This could be done by using new efficient finite-difference 
methods or by using the mixing theorem. The latter means that conditional on the 
whole path of the instantaneous variance from vq to vt the option price G{ then is given 
by the Black-Scholes formula with the efficient spot price S e and the efficient volatility 
a e (see Romano and Touzi (1997)). This method is very fast because its complexity is 



equivalent to one-dimensional MC. 
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